A Simple Method for Computing the Non-Linear Mass 
Correlation Function with Implications for Stable Clustering 



Robert R. Caldwell^'^, Roman Juszkiewicz^ ^, 
Paul J. Steinhardt^'^, and Frangois R. Bouchet^'^ 

ABSTRACT 

We propose a simple and accurate method for computing analytically the mass 
correlation function for cold dark matter and scale-free models that fits N-body 
simulations over a range that extends from the linear to the strongly non-linear 
regime. The method, based on the dynamical evolution of the pair conservation 
equation, relies on a universal relation between the pair-wise velocity and the 
smoothed correlation function valid for high and low density models, as derived 
empirically from N-body simulations. An intriguing alternative relation, based 
on the stable-clustering hypothesis, predicts a power-law behavior of the mass 
correlation function that disagrees with N-body simulations but conforms well 
to the observed galaxy correlation function if negligible bias is assumed. The 
method is a useful tool for rapidly exploring a wide span of models and, at the 
same time, raises new questions about large scale structure formation. 
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Understanding the origin and evolution of the clustering pattern of galaxies is one of 
the most important goals of cosmology. Until now, this problem has been investigated using 
a four- fold path: (1) perturbation theory (for a review of recent advances, see [Juszkiewicz fc 



[Bouchet 1996| and references therein); (2) a kinetic description, adapted from the BBGKY 
hierarchy, used in plasma physics (Peebles 19801 , §IV); (3) N-body simulations (e.g., [Jenkins 



\et al. 19"9^ , hereafter VIRGO); (4) and semi-analytic fits to N-body results, based on the 
so-called universal scaling hypothesis (see [Hamilton et al. 1991| , [Jain et al. 1995[ , [Peacock 
fc Dodds 199^ , [Ma 1998| 
plementary. 



The advantages and limitations of these methods are often com- 
For example, applying perturbation theory often leads to analytic results for 
a wide class of models while the N-body simulations allow a study of only one model at a 
time. On the other hand, perturbation theory works only in the weakly non-linear regime 
while N-body experiments describe the fully non-linear dynamics, albeit over a limited dy- 
namical range. The subject of the present study is an analytic ansatz for the evolution of 
the two-point correlation function of density fluctuations spanning the linear and non-linear 
regime, which builds on all four methods described above. 

Our approach is based on the pair conservation equation, which relates the mean (pair- 
weighted) relative velocity of a pair of particles to the time evolution of the correlation 
function in a self-gravitating gas: 

3[l + ^{x,a)] da ~ Hr ' ^ ' 



where a(t) is the expansion factor with a = 1 at present, r = ax is the proper separation 
and H{a) is the Hubble parameter (see pavis fc Peebles 1977[ , [Peebles 1980|) . Here ^{x,a) 
represents the two-point correlation function averaged over a ball of comoving radius x: 



3 

x-^ Jo 



(2) 



The approximate solution of (|I]) is known in the large separation limit, where |^| -C 1 (linear 
regime); the stable clustering hypothesis is often invoked to describe the small separation 
limit, where ^ ^ 1 (non-linear regime). Hence, equation (|I]) is "a guide to speculation on the 
behavior of the correlation function" ( [Peebles 1980| , p. 268) since an assumed implies a 
function ^ that should agree with the weak and strong fleld limits, and interpolates between 
these limits in a reasonable way. 

An approximate universal relation between the pair-wise velocity and the smoothed 
correlation function has been conjectured by [Hamilton et al. 1991| on the basis of N-body 



simulation results and further explored in [Nityananda fc Padmanabhan 199^ and [Padman- 
abhan fc Engineer 1998[ . In the past, the relation has been used in attempts to derive a 
general functional that converts directly from a linear to a non-linear mass correlation. In 
this paper, we present a simple extension of the relation that applies to both high- and 
low-density models, but take a different approach to obtaining the non-linear correlation 
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function. Namely, we use the universal relation to close (|l|); we then evolve the resulting 
partial differential equation to compute the non-linear correlation function. This turns out 
to be a fast and surprisingly accurate method that matches N-body results for a wide variety 
of cold dark matter (CDM) models. As a stand-alone computer program, the algorithm can 
be adopted on a programmable calculator; or, it can easily be incorporated in more sophis- 
ticated programs that predict other cosmological properties, such as CMBFAST ( [Seljak fc| 
Zaldarriaga 1996|) . 



Figure |1| clearly illustrates the nearly model-independent relationship between the pair- 
wise velocity and the smoothed correlation function, observed in N-body simulations for a 
wide range of perturbation spectra. We define this relation as 

nf^]^-^- (3) 



Compared to [Hamilton et at. 1991| , a novel feature is plotting the relation in terms of 



where / = d\n.D / d\n.a is the standard linear density growth factor, rather than ^ alone — 
a difference that is essential for extending the relation to low-density models. The evolved, 
non-linear clustering of scale-free spectra with = — 1, —2 as well as the CDM family of 
models produces a very similar relation between —v^2/Hr and ^. An excellent fit to the 
functional relation V[fE,{x, a)] in Figure |l], based on the n = — 1 curve, is given by 

{\x X < 0.15 

0.7x exp(-0.31x°-^^) 0.15 < x < 20 (4) 
3.3x-°-i^ x>20 

valid for x ~ 10^. In this paper, we use this fitting formula, designed for the n = — 1 curve, 
as the expression for V[x\ in (|1]) to be applied to all models. We find this to be sufficient to 
reproduce N-body results to within 10% accuracy over the range of models and scales shown 
in Figure ^ which extends deep into the non-linear regime. To push to lower density models 
and improve the accuracy further, it would be simple to modify the algorithm for V[x] to 
include the ^-dependent rise near x ~ 20. 

Starting with the linear correlation function, and armed only with information about 
the background cosmological evolution, we propose to dynamically obtain the non-linear 
as a function of separation and time. Here then is our idealized procedure in three steps. 

1. Reformulate: We first rewrite the partial differential equation as 

- ^'^^vm (5) 

a In a 4 
where V[x\ is given in equation (|). 

2. Initialize: The initial conditions are set by the linear correlation function at a red shift 
z = — 1 + 1/aj such that ^{x, ai) is less than unity for all x of interest. Our procedure assumes 
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that only the amplitude, and not the shape of the linear correlation function has changed over 
this interval, as occurs for cold dark matter models with a primordial spectrum of adiabatic 
density perturbations. Hence, this procedure will not apply to cosmological models with a 
late-time decaying neutrino, but will apply to hot dark matter models wherein the shape of 
the linear power spectrum is set by red shift z ~ 100. 

3. Evolve: We numerically solve the partial differential equation and evolve ^. At each 
step in the evolution, we use ^ = ^ x (1 — 7/3) with 7 = — dln^/dlnx to determine the 
correlation function. The value of / is updated at each step in a as appropriate for the 
cosmology. 

The remarkable results are shown in Figure |^. Here we see that our simple procedure 
gives excellent agreement with N-body simulations. Based on the span of behavior in the 
cosmic time evolution and the shape of correlation function, we expect this procedure should 
be valid for a wide range of cosmological models, including quintessence ([Caldwell, Dave,"fc 
Steinhardt 1998| ) and models with tilted spectra. 



Figure |] demonstrates that we have obtained a simple and powerful new tool for rapidly 
and accurately obtaining the non-linear power spectrum for a wide range of models. How- 
ever, the physical origin of the nearly model-independent relation is not understood 
in detail. In the linear regime, perturbation theory predicts —VyijHr = (2/3)/^. In the 
non-linear regime, Padmanabhan et al. (1996) have suggested that insight may be obtained 
by comparison to the case of the gravitational collapse of a spherical top hat mass distri- 
bution. Using their solution (eqs. (16-19) in their paper), we find that —v^^jHr is linearly 
proportional to times a slowly decreasing function of ^ for a surprisingly wide range 
of ^ ^ 1, including ^ = 10, the turnover point in Figure |l[ In particular, for ^ = 10, 
— fia/i/r = 1.77 /(fi), which is similar to — fia/i^r ^ 2f{Q) in Figure |1[ 

In the strongly non-linear regime, ^ ^ 10, Figure shows a visible difference in the 
shape of ^[/^] between the high density and low density models. This may be due to the 
suppression of linear growth, which occurs at late times in low density models and leads 
to the enhanced clustering on small scales relative to large scales. However, this has a 
negligible effect on the computed non-linear correlation function. For example, using the 
curve for ACDM shown in Figure |^ as the basis of V in our procedure, we find the amplitude 
of the non-linear correlation function differs by only 10% at r ~ 0.1 Mpc/h. For the models 
shown, this accuracy is comparable to what is obtained by Hamilton et al. 1991| and Peacock 



fc Dodds 1996|. The advantage here is that our method can be immediately applied to new 



types of CDM models {e.g., quintessence cosmologies) without having to run new N-body 
simulations to fix fitting parameters. 

An important issue raised by our ansatz is the validity of the stable clustering hypothesis. 
The stable clustering regime corresponds to the limit where particle pairs detach from the 
Hubble flow and —Vi2/Hr — > 1. Figure |l| shows that —Vii/Hr first overshoots unity by a 
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factor of two and then rebounds towards unity. However, it is not clear whether it converges 
to unity at ,^ ^ 1000 or possibly oscillates if the simulations are extended to higher values 

ofe 

It is interesting to compare the predictions of our ansatz if the relation between —v^2/Hr 
and ^ is modified to enforce more rapid convergence to stable clustering. A ready exam- 
ple is an alternative ansatz based on the pair conservation equation recently proposed by 
[Juszkiewicz et al. 1999| (hereafter JSD). Their ansatz for f^a, based on an interpolation 
between the behavior predicted by perturbation theory in the weakly non-linear regime and 
stable clustering in the strongly non-linear regime is given by 

fi2(a;, a) := -^Hrf^{x,a) 1 + a^{x,a) (6) 

where ^ = ^/{l + C,) and a is a function which controls the strength of the non-linear 
feedback. Here we use a = 1.8 — I.I7, based on perturbation theory (see [Scoccimarro fc 



b'rieman 199(j| ), where 7 is the slope of the correlation function at ^ = 1. The key point, as 



shown in Figure y, is that the pair-wise velocity rapidly approaches the stable clustering limit 
by ^ ~ 10, and remains there to within ~ 20% on smaller scales in the more strongly non- 
linear regime. This means particle pairs separated by ^ 1 Mpc/h have the rough behavior 
of virialized objects, such as clusters and galaxies. 

The correlation function ^ obtained by closing the pair conservation equation with (|^), 
as shown in Figure ^, displays a power-law behavior in the non-linear regime with index 
~ —1.7, in disagreement with N-body simulations of CDM but curiously similar to the 
galaxy correlation function observed in the APM survey ([Maddox et al. 1996| ). This result is 
not unique to the JSD ansatz; substituting any shape similar to that shown in the top panel 
of Figure |^ for into our ansatz would produce a similar effect on the mass correlation 

function. In the past, the disagreement between the dark matter power spectrum observed 
in simulations, which shows no evidence of power-law behavior, and the simple power-law 
observed in the galaxy correlation function has been attributed to a scale-dependent bias. 
The conventional picture is that the bias function, b{r), is just so as to cause a non-power- 
law behavior in the dark matter to be translated into a power-law behavior of the galaxy 
correlation function, C,g{r) = ^(r)6^(r). 

Our present findings concerning the dependence of the mass correlation function on 
the model independent —v^2/Hr vs. relation suggests a radical possibility. Perhaps the 
bias factor is completely negligible and, instead, CDM models are missing some important 
physical feature {e.g. mechanics of galaxy formation, or some property of the dark matter) 
which causes rapid convergence to stable clustering in the non-linear regime {—V12/ Hr — > 1 
without substantial overshoot) and to power-law behavior of the correlation function. The 
ansatz illustrated in Figure ^, which enforces stable clustering by fiat, may be implicitly 
describing a modified CDM model of galaxy formation which incorporates the new physical 
feature. The difference between N-body simulations and observations, whether due to bias 
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in the conventional picture or to something more radical, such as a modification to CDM, 
can perhaps be determined empirically by studying red shift distortion on the scales where 
N-body simulation suggests overshoot in —Vy^jHr and the ansatz of Figure ^ does not. 

In sum, our studies have produced a simple recipe for computing the non-linear power 
spectrum for a wide range of models. (Upon publication of the paper, we will make a program 
available at feynman.princeton.edu/~steinh.) A key feature of the method is the universal 
function relating the pair velocity to the mass correlation function which does not converge 
rapidly to the stable clustering limit, but, rather, overshoots by a factor of two. This feature 
is responsible for the fact that the mass correlation function does not approach a power-law. 
Our studies have also raised several interesting issues in structure formation: why is —v^^jHr 
vs. so similar for a wide range of models? can the universal relation be derived from 
theory? for what range of models is the relation model-independent? how does the universal 
relation ultimately approach the stable clustering limit at small scales (if it does at all)? 
and, does the success of a universal relation based on the stable clustering hypothesis, as in 
Figure 0, suggest a viable, alternative explanation for the power-law behavior of the galaxy 
correlation function? 
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10"' 10° 10' _ lO" 10^ 10^ 

Fig. 1. — The pair-wise velocity in units of the Hubble velocity, — Wia/i^r is shown as a 
function of the smoothed correlation function times the linear density growth factor, / ^ as 
determined by N-body simulations (scale free and SCDM from |Jain 1997] , Figures 3 & 8; A 
and OCDM from VIRGO). Not only for these examples, but for a wider variety of initial 
conditions and at different times, the pair-wise velocity displays nearly the same shape in ^, 
leading us to a nearly universal relation, V[f^] = —v^2/Hr. 
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Fig. 2. — The non-linear correlation function ^ as a function of separation is shown for various 
cases. The non-linear ^ obtained from our method (solid lines) is in excellent agreement with 
the N-body results (circles). The dashed line is the linear ^. The cosmological models are 
from the VIRGO simulation, whereas the n = —1 scale-free model is from | Jain 1997| , Figure 
1. For n = —1, X = xnl at ^ = 1. 
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Fig. 3. — The pair-wise velocity and non-linear correlation function are shown for a case (solid 
curves) in which we have forced a rapid convergence to stable clustering, —v^2/Hr = 1. 
For the purposes of illustration, we have used the JSD ansatz described by equation (|^), 
although substituting any similar shape for would produce a similar effect on the 

mass correlation function. The rapid approach to stable clustering and the absence of any 
significant overshoot, as seen in the solid curve on the left, results in the power law behavior 
^ oc r~^'^ on the right. The circles and squares in each panel are the N-body results for 
SCDM. 



